library(Seurat)
library(spacexr)
library(tidyverse)
library(spatstat)
library(nlme)
library(lgcp)
library(CAinterprTools)
library(sf)
library(concaveman)
library(ggplot2)

loading dataset

load(file= "/Users/cuihening/Desktop/shikun/seurat/robjs/chicken_visium.4.Robj")
TM <- FetchData(object=chicken_visium,vars=c('TMSB4X','ident', 'orig.ident')) 
tmsb4x=cbind(cell = rownames(TM), TM) %>% 
  dplyr::rename(time = orig.ident) 
tmsb4x$is2=as.factor(ifelse(tmsb4x$ident==2,1,0))

The function to get coordinate

getcoor=function(slice, time){
coords <- GetTissueCoordinates(chicken_visium, image= slice)
colnames(coords) <- c("x", "y")
coords <- cbind(cell = rownames(coords),coords)
coords$time = time
return(coords)
}
coords_d4 = getcoor('slice1','D4')
coords_d10 = getcoor('slice1_D10-C1','D10')
coords_d7 = getcoor('slice1_D7-B1','D7')
coords_d14 = getcoor('slice1_D14-D1','D14')

Function for constuct ppp

buildppp=function(dataset,markers){
X <- dataset$x
Y <- dataset$y
xrange <- range(X, na.rm=T)
print(xrange)
yrange <- range(Y, na.rm=T)
print(yrange)
ppppro =  ppp(X,Y,xrange,yrange,marks=markers)
return(ppppro)
}

EDA

D4

D4 = tmsb4x %>% 
  filter(time=='D4') %>% 
  merge(coords_d4, by="cell")  %>% 
  filter(x >= 385 &  y>=267 )  %>% 
  rename(Y=x, X=y) %>% 
  rename(x=X, y=Y) %>% 
  mutate(x=x-138, y=y-208) %>% 
  mutate(time=4)
summary(D4)
##      cell               TMSB4X            ident       time.x          is2    
##  Length:141         Min.   :0.04300   6      :66   Length:141         0:119  
##  Class :character   1st Qu.:0.05200   4      :27   Class :character   1: 22  
##  Mode  :character   Median :0.05554   2      :22   Mode  :character          
##                     Mean   :0.05688   8      :21                             
##                     3rd Qu.:0.05953   5      : 5                             
##                     Max.   :0.07820   0      : 0                             
##                                       (Other): 0                             
##        y               x            time.y               time  
##  Min.   :182.6   Min.   :142.5   Length:141         Min.   :4  
##  1st Qu.:202.2   1st Qu.:180.0   Class :character   1st Qu.:4  
##  Median :221.8   Median :206.3   Mode  :character   Median :4  
##  Mean   :222.8   Mean   :201.6                      Mean   :4  
##  3rd Qu.:241.3   3rd Qu.:225.1                      3rd Qu.:4  
##  Max.   :274.0   Max.   :247.5                      Max.   :4  
## 
xym= D4 %>% select(c("x","y")) 
pts <- st_as_sf(xym, coords=c('x','y'))
conc <- concaveman(pts)
ggplot() +
  geom_sf(data = conc, fill = NA)

d4_owin= as.owin(as_Spatial(conc))
plot(d4_owin)

D4_samppp =  ppp(D4$x,D4$y, c(0, 450),c(0,550), marks=D4$is2)
D4_samppp =  ppp(D4$x,D4$y, window = d4_owin, marks=D4$is2)
## Warning: 19 points were rejected as lying outside the specified window
D4_ty2=D4 %>% 
  filter(is2==1)
D4_ty2p=ppp(D4_ty2$x,D4_ty2$y, c(0, 450),c(0,550))
D4_ty2_owin=ppp(D4_ty2$x,D4_ty2$y, window = d4_owin)
## Warning: 4 points were rejected as lying outside the specified window
par(mfrow=c(1,2))
plot(D4_samppp, pch = 20,show.window=FALSE,cols=c("blue", "green"), legend=FALSE, main="D4 point",markscale = 0.1)
## Warning in plot.ppp(D4_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 19 illegal points also plotted
plot(density(D4_ty2_owin), main="D4 density")

pixel=pixellate.ppp(D4_ty2_owin,fractional=TRUE)
plot(pixel)

plot(as.im(density(D4_ty2_owin), fractional=TRUE))

par(mfrow=c(1,2))
plot(D4_samppp, pch = 20,show.window=FALSE,cols=c("blue", "green"), legend=FALSE, main="D4 point",markscale = 0.1)
## Warning in plot.ppp(D4_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 19 illegal points also plotted
plot(density(D4_ty2p), main="D4 density")

D4_fit =  kppm(D4_ty2p ~1, "LGCP")
D4_fit
## Stationary Cox point process model
## Fitted to point pattern dataset 'D4_ty2p'
## Fitted by minimum contrast
##  Summary statistic: K-function
## 
## Uniform intensity:   8.888889e-05
## 
## Cox model: log-Gaussian Cox process
##  Covariance model: exponential
## Fitted covariance parameters:
##       var     scale 
##  5.041248 62.154818 
## Fitted mean of log of random intensity: -11.84875
D4_fit =  kppm(D4_ty2_owin ~1, "LGCP")
D4_fit
## Stationary Cox point process model
## Fitted to point pattern dataset 'D4_ty2_owin'
## Fitted by minimum contrast
##  Summary statistic: K-function
## 
## Uniform intensity:   0.002988242
## 
## Cox model: log-Gaussian Cox process
##  Covariance model: exponential
## Fitted covariance parameters:
##          var        scale 
## 1.036636e+00 3.273900e+07 
## Fitted mean of log of random intensity: -6.331388
plot(predict.dppm(D4_fit))

D7

D7 = tmsb4x %>% 
  filter(time=='D7') %>% 
  merge(coords_d7, by="cell")  %>% 
  filter(x >= 255 &  y>=267 )  %>% 
  rename(Y=x, X=y) %>% 
  rename(x=X, y=Y) %>% 
  mutate(x=x-150, y=y-152) %>% 
  mutate(time=7)
summary(D7)
##      cell               TMSB4X            ident        time.x          is2    
##  Length:505         Min.   :0.00000   1      :166   Length:505         0:438  
##  Class :character   1st Qu.:0.05050   5      : 83   Class :character   1: 67  
##  Mode  :character   Median :0.05403   2      : 67   Mode  :character          
##                     Mean   :0.05413   4      : 53                             
##                     3rd Qu.:0.05735   8      : 48                             
##                     Max.   :0.07924   7      : 35                             
##                                       (Other): 53                             
##        y               x            time.y               time  
##  Min.   :107.8   Min.   :131.1   Length:505         Min.   :7  
##  1st Qu.:179.5   1st Qu.:172.2   Class :character   1st Qu.:7  
##  Median :225.1   Median :198.6   Mode  :character   Median :7  
##  Mean   :220.9   Mean   :200.3                      Mean   :7  
##  3rd Qu.:264.3   3rd Qu.:228.6                      3rd Qu.:7  
##  Max.   :323.1   Max.   :281.1                      Max.   :7  
## 
xym= D7 %>% select(c("x","y")) 
pts <- st_as_sf(xym, coords=c('x','y'))
conc <- concaveman(pts)
ggplot() +
  geom_sf(data = conc, fill = NA)

d7_owin= as.owin(as_Spatial(conc))
plot(d7_owin)

D7_samppp =  ppp(D7$x,D7$y, window = d7_owin, marks=D7$is2)
## Warning: 34 points were rejected as lying outside the specified window
D7_ty2=D7 %>% 
  filter(is2==1)
D7_ty2p=ppp(D7_ty2$x,D7_ty2$y, window = d7_owin)
## Warning: 3 points were rejected as lying outside the specified window
par(mfrow=c(1,2))
plot(D7_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D7 point")
## Warning in plot.ppp(D7_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 34 illegal points also plotted
plot(density(D7_ty2p), main="D7 density")

D7_fit =  kppm(D7_ty2p ~1, "LGCP")
D7_fit
## Stationary Cox point process model
## Fitted to point pattern dataset 'D7_ty2p'
## Fitted by minimum contrast
##  Summary statistic: K-function
## 
## Uniform intensity:   0.002783736
## 
## Cox model: log-Gaussian Cox process
##  Covariance model: exponential
## Fitted covariance parameters:
##          var        scale 
## 8.909741e-01 1.332978e+08 
## Fitted mean of log of random intensity: -6.329448

D10

D10 = tmsb4x %>% 
  filter(time=='D10') %>% 
  merge(coords_d10, by="cell")  %>% 
  filter(x >= 193 &  y>=320)  %>% 
  rename(Y=x, X=y) %>% 
  rename(x=X, y=Y) %>% 
  mutate(x=x-217, y=y-128) %>% 
  mutate(time=10)
summary(D10)
##      cell               TMSB4X            ident        time.x          is2    
##  Length:888         Min.   :0.03259   3      :163   Length:888         0:760  
##  Class :character   1st Qu.:0.05015   2      :128   Class :character   1:128  
##  Mode  :character   Median :0.05414   1      :123   Mode  :character          
##                     Mean   :0.05418   4      :118                             
##                     3rd Qu.:0.05806   0      :109                             
##                     Max.   :0.07277   7      :108                             
##                                       (Other):139                             
##        y                x            time.y               time   
##  Min.   : 72.12   Min.   :103.4   Length:888         Min.   :10  
##  1st Qu.:169.94   1st Qu.:156.0   Class :character   1st Qu.:10  
##  Median :215.79   Median :197.4   Mode  :character   Median :10  
##  Mean   :220.37   Mean   :200.8                      Mean   :10  
##  3rd Qu.:267.97   3rd Qu.:238.5                      3rd Qu.:10  
##  Max.   :378.87   Max.   :321.0                      Max.   :10  
## 
xym= D10 %>% select(c("x","y")) 
pts <- st_as_sf(xym, coords=c('x','y'))
conc <- concaveman(pts)
ggplot() +
  geom_sf(data = conc, fill = NA)

d10_owin= as.owin(as_Spatial(conc))
plot(d10_owin)

D10_samppp =ppp(D10$x,D10$y, window = d10_owin, marks=D10$is2)
## Warning: 46 points were rejected as lying outside the specified window
D10_ty2=D10 %>% 
  filter(is2==1)
D10_ty2p=ppp(D10_ty2$x,D10_ty2$y, window = d10_owin)
## Warning: 2 points were rejected as lying outside the specified window
par(mfrow=c(1,2))
plot(D10_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D10 point")
## Warning in plot.ppp(D10_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 46 illegal points also plotted
plot(density(D10_ty2p), main="D10 density")

D10_fit =  kppm(D10_ty2p ~1, "LGCP")
D10_fit
## Stationary Cox point process model
## Fitted to point pattern dataset 'D10_ty2p'
## Fitted by minimum contrast
##  Summary statistic: K-function
## 
## Uniform intensity:   0.002817653
## 
## Cox model: log-Gaussian Cox process
##  Covariance model: exponential
## Fitted covariance parameters:
##          var        scale 
## 1.000520e+00 1.545034e+06 
## Fitted mean of log of random intensity: -6.372111

D14

D14 = tmsb4x %>% 
  filter(time=='D14') %>% 
  merge(coords_d14, by="cell")  %>% 
  rename(Y=x, X=y) %>% 
  rename(x=X, y=Y) %>% 
  mutate(x=x-144, y=y-78) %>% 
  mutate(time=14)
summary(D14)
##      cell               TMSB4X            ident        time.x          is2     
##  Length:1967        Min.   :0.03193   0      :854   Length:1967        0:1643  
##  Class :character   1st Qu.:0.05184   2      :324   Class :character   1: 324  
##  Mode  :character   Median :0.05530   3      :175   Mode  :character           
##                     Mean   :0.05509   4      :159                              
##                     3rd Qu.:0.05871   8      :138                              
##                     Max.   :0.07352   7      :128                              
##                                       (Other):189                              
##        y                 x             time.y               time   
##  Min.   :  9.428   Min.   : 58.58   Length:1967        Min.   :14  
##  1st Qu.:139.719   1st Qu.:137.43   Class :character   1st Qu.:14  
##  Median :224.323   Median :200.63   Mode  :character   Median :14  
##  Mean   :220.551   Mean   :200.92                      Mean   :14  
##  3rd Qu.:302.755   3rd Qu.:261.00                      3rd Qu.:14  
##  Max.   :400.476   Max.   :361.98                      Max.   :14  
## 
xym= D14 %>% select(c("x","y")) 
pts <- st_as_sf(xym, coords=c('x','y'))
conc <- concaveman(pts)
ggplot() +
  geom_sf(data = conc, fill = NA)

d14_owin= as.owin(as_Spatial(conc))
plot(d14_owin)

D14_samppp =  ppp(D14$x,D14$y, window = d14_owin, marks=D14$is2)
## Warning: 64 points were rejected as lying outside the specified window
D14_ty2=D14 %>% 
  filter(is2==1)
D14_ty2p=ppp(D14_ty2$x,D14_ty2$y, window = d14_owin)
## Warning: 11 points were rejected as lying outside the specified window
par(mfrow=c(1,2))
plot(D14_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D14 point")
## Warning in plot.ppp(D14_samppp, pch = 20, show.window = FALSE, cols = c("blue",
## : 64 illegal points also plotted
plot(density(D14_ty2p), main="D14 density")

D14_fit =  kppm(D14_ty2p ~1, "LGCP")
D14_fit
## Stationary Cox point process model
## Fitted to point pattern dataset 'D14_ty2p'
## Fitted by minimum contrast
##  Summary statistic: K-function
## 
## Uniform intensity:   0.00337893
## 
## Cox model: log-Gaussian Cox process
##  Covariance model: exponential
## Fitted covariance parameters:
##         var       scale 
##    1.216153 4084.813883 
## Fitted mean of log of random intensity: -6.298273

Spatstat replicate model

Create the list of all point

time_ty2=solist(D4_ty2_owin,D7_ty2p, D10_ty2p, D14_ty2p)
time_ty2
## List of point patterns
## 
## Component 1:
## Planar point pattern: 18 points
## window: polygonal boundary
## enclosing rectangle: [142.5205, 247.5368] x [182.562, 274.0155] units
## *** 4 illegal points stored in attr(,"rejects") ***
## 
## Component 2:
## Planar point pattern: 64 points
## window: polygonal boundary
## enclosing rectangle: [131.0609, 281.1149] x [107.7713, 323.0721] units
## *** 3 illegal points stored in attr(,"rejects") ***
## 
## Component 3:
## Planar point pattern: 126 points
## window: polygonal boundary
## enclosing rectangle: [103.3981, 320.9683] x [72.1171, 378.8714] units
## *** 2 illegal points stored in attr(,"rejects") ***
## 
## Component 4:
## Planar point pattern: 313 points
## window: polygonal boundary
## enclosing rectangle: [58.5757, 361.9798] x [9.428, 400.4762] units
## *** 11 illegal points stored in attr(,"rejects") ***
sapply(time_ty2, npoints)
## [1]  18  64 126 313
plot.solist(time_ty2, main="The ty2 cell", main.panel=letters[1:4], legend=FALSE, equal.scales = TRUE, halign=TRUE, valign=TRUE)

Turn list into hyperframe

hyoertime <- hyperframe(cell=time_ty2)
hyoertime$time=c(4,7,10,14)
plot(hyoertime, quote(plot(density(cell), main=time)))

try=mppm(cell~time, hyoertime)

test the model

cdf.test.mppm(try, "y")
## 
##  Spatial Kolmogorov-Smirnov test
## 
## data:  predicted cdf of covariate '"y"' evaluated at data points of 'try'
## D = 0.59049, p-value < 2.2e-16
## alternative hypothesis: two-sided
plot.mppm(try, main="The ty2 cell", cif=TRUE,
                how="image")
## Warning: glm.fit: algorithm did not converge

lurking(try, covariate = expression(time))
## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge

sub_time=subfits(try)
## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: algorithm did not converge
fitted_intensities <- predict.mppm(try)
fitted_intensities2 <- predict.mppm(try,locations = D4_ty2_owin)
fitted_intensities2
## Hyperframe:
##   trend  cif
## 1  (im) (im)
## 2  (im) (im)
## 3  (im) (im)
## 4  (im) (im)
K_functions4 <-  Kinhom(D4_ty2_owin, ratio= TRUE)
K_functions7 <-  Kinhom(D7_ty2p, ratio= TRUE)
K_functions10 <-  Kinhom(D10_ty2p, ratio= TRUE)
K_functions14 <-  Kinhom(D14_ty2p, ratio= TRUE)
Kcomb=pool(K_functions4,K_functions7,K_functions10,K_functions14)
try=lgcp.estK(Kcomb)
try
## Minimum contrast fit (object of class "minconfit")
## Model: log-Gaussian Cox process
##   Covariance model: exponential
## Fitted by matching theoretical K[inhom] function to Kcomb
## 
## Internal parameters fitted by minimum contrast ($par):
##       sigma2        alpha 
## 1.2216365337 0.0001424789 
## 
## Fitted covariance parameters:
##          var        scale 
## 1.2216365337 0.0001424789 
## 
## Converged successfully after 57 function evaluations
## 
## Starting values of parameters:
## sigma2  alpha 
##      1      1 
## Domain of integration: [ 0 , 22.86 ]
## Exponents: p= 2, q= 0.25

Cox Spatial Temporal Model

prepare data

temp= rbind(D4, D7, D10, D14) %>% 
  filter(ident==2) %>% 
  select(x,y,time, is2)
temppx=ppx(temp)
plot(temppx)

library("stpp")
## Warning: package 'stpp' was built under R version 4.1.2
## Loading required package: rpanel
## Warning: package 'rpanel' was built under R version 4.1.2
## Loading required package: tcltk
## Package `rpanel', version 1.1-5: type help(rpanel) for summary information
## 
## Attaching package: 'rpanel'
## The following object is masked from 'package:tidyr':
## 
##     population
## Loading required package: splancs
## Warning: package 'splancs' was built under R version 4.1.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.1.2
## 
## Spatial Point Pattern Analysis Code in S-Plus
##  
##  Version 2 - Spatial and Space-Time analysis
## 
## Attaching package: 'splancs'
## The following object is masked from 'package:dplyr':
## 
##     tribble
## The following object is masked from 'package:tidyr':
## 
##     tribble
## The following object is masked from 'package:tibble':
## 
##     tribble
tlim=c(0,14)
xyt <- stppp(list(data = temp, tlim = tlim, window = owin(c(0,450),c(0,550))))
xyt <- integerise(xyt)
xyt
## Space-time point pattern
## Planar point pattern: 541 points
## window: rectangle = [0, 450] x [0, 550] units
##    Time Window : [ 0 , 14 ]
den <- lambdaEst(xyt, axes = TRUE)
plot(den)
sar <- spatialAtRisk(den)
sar
## SpatialAtRisk object
##    X range : [1.7578125,448.2421875]
##    Y range : [2.1484375,547.8515625]
##    dim     : 128 x 128
plot(sar)

mut <- muEst(xyt)
mut
## temporalAtRisk object
## function (t) 
## {
##     if (length(t) > 1) {
##         stop("Function only works for scalar numeric objects, ie a vector of length 1.")
##     }
##     if (!any(as.integer(t) == tvec)) {
##         return(NA)
##     }
##     return(obj[which(as.integer(t) == tvec)] * scale)
## }
## <bytecode: 0x7fa201af6240>
## <environment: 0x7fa201af7468>
## attr(,"tlim")
## [1]  0 14
## attr(,"class")
## [1] "temporalAtRisk" "function"      
##    Time Window : [ 0 , 14 ]
plot(mut)

gin <- ginhomAverage(xyt, spatial.intensity = sar, temporal.intensity = mut)
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |======================================================================| 100%
## Returning an average of 4 curves
plot(gin)

kin <- KinhomAverage(xyt, spatial.intensity = sar, temporal.intensity = mut)
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |======================================================================| 100%
## Returning an average of 4 curves